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EXECUTIVE  SUMMARY 


Sea  clutter  mitigation  is  an  essential  feature  of  airborne  radars  operating  in  a  maritime  environment 
because  the  radar  return  from  the  sea  surface  can  easily  obscure  targets  of  interest,  therefore  being  able  to 
accurately  model  the  sea  clutter  returns  is  an  essential  step  in  the  design  process. 

Much  work  has  been  done  in  the  area  of  sea  clutter  modeling,  including  measurement  campaigns,  sta¬ 
tistical  characterization,  and  electromagnetic  (EM)  scattering  computations.  However,  these  models  arc 
usually  purpose-built  to  solve  a  specific  problem  and  arc  rarely  general  in  character.  A  more  comprehensive 
statistical  characterization  is  also  desired  that  does  not  impose  a  large  computational  burden  like  an  EM 
scattering  computation  would. 

In  this  work  we  present  a  general  framework  for  modeling  sea  clutter  returns  from  a  signal  processing 
perspective  by  modeling  the  sea  clutter  scattering  channel  as  a  random  linear  time- varying  system.  Linear 
time- varying  system  theory  has  been  studied  rigorously  since  the  1960's,  and  by  casting  the  clutter  scattering 
problem  in  this  form,  we  can  apply  standard  analytical  tools  immediately  to  gain  new  insight  into  the  time- 
frequency  characteristics  of  clutter  returns. 

The  mathematical  formulation  of  sea  clutter  in  this  work  is  hoped  to  provide  a  unifying  structure  that 
will  enhance  the  future  development  of  clutter  mitigation  techniques,  particularly  in  the  design  of  radar 
waveforms  and  signal  processing  techniques. 


E-l 


APPLICATION  OF  STATISTICAL  LINEAR  TIME-VARYING  SYSTEM  THEORY  TO 
MODELING  OF  HIGH  GRAZING  ANGLE  SEA  CLUTTER 


1.  INTRODUCTION 
1.1  Proposed  Approach 

Another  approach  to  modeling  clutter  which  has  not  received  much  attention  in  the  literature  is  to  use 
tools  from  linear  time-varying  system  theory.  A  linear  time-varying  (LTV)  system  is  one  whose  properties 
change  with  time.  In  the  case  of  a  linear  time-invariant  (LTI)  system  the  output  signal  y(t)  can  be  expressed 
using  the  convolution  of  the  input  signal  x(t)  with  the  impulse  response  h(t): 

y(t )  =  J  h{r)x(t  -  t )dr.  (1) 

In  this  case,  y(t)  can  be  seen  as  the  superposition  of  copies  of  x(t),  each  delayed  by  a  lag  r  and  scaled  by  a 
coefficient  h(r).  In  an  LTV  system,  the  impulse  response  depends  on  both  t  and  r  [1]  and  thus  the  scaling 
coefficient  now  changes  with  time  and  becomes  h(r,  t): 

y(t)  =  J  h(r,  t)x(t  -  t)(It.  (2) 

Due  to  the  fact  that  the  LTV  approach  is  very  prevalent  in  wireless  communication  theory  in  modeling 
multipath  fading  channels,  it  is  common  to  refer  to  the  system  represented  by  h  as  the  “channel”.  Changes 
in  h  with  t  characterize  changes  in  the  channel  over  time,  e.g.  if  a  scatterer  is  moving,  which  can  introduce 
Doppler  frequency  shifts  in  the  received  signal.  (In  this  work  we  assume  the  velocities  involved  are  small 
relative  to  the  speed  of  light  and  thus  that  the  Doppler  effect  affects  only  the  carrier  frequency  and  not  the 
time-scaling  of  the  signal  envelope.) 

Because  of  the  multitude  of  small  variations  on  the  sea  surface  and  the  impracticality  of  modeling  the 
exact  sea  surface  at  a  specified  point  on  the  earth,  we  will  model  the  radar  reflection  off  the  sea  surface  as  a 
random  linear  time-varying  system.  The  fundamental  work  on  random  LTV  system  theory  was  performed 
by  Bello  [2]  where  he  introduces  the  concept  of  a  wide-sense  stationary  uncorrelated  scattering  (WSSUS) 
channel. 

LTV  system  theory  has  been  primarily  used  by  wireless  communications  researchers  [1,  3,  4]  for  de¬ 
signing  adaptive  radios  that  adjust  their  modulations  based  on  changing  channel  parameters,  such  as  the 
channel  coherence  time  and  coherence  bandwidth,  which  describe  properties  of  h(r,t).  However,  one  of 
the  initial  uses  of  random  LTV  system  theory  was  in  describing  distributed  radar  returns  from  a  large  object 
(such  as  a  planet)  in  radio  astronomy  as  a  function  of  delay  and  Doppler  shift  using  a  so-called  “scattering 
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function”  [5,6].  The  use  of  scattering  functions  for  modeling  radar  returns  arc  also  described  by  Van  Trees 
in  his  seminal  book  series  [7],  by  Kay  in  at  least  one  paper  [8],  and  are  mentioned  by  a  small  number  of 
other  authors  (e.g.  [9-15]),  but  in  general  arc  not  commonly  referred  to  in  the  radar  literature  or  used  widely 
among  practicing  radar  engineers. 

In  this  work  we  propose  the  use  of  LTV  system  theory  for  modeling  sea  clutter  returns  so  that  these  tools 
from  decades  of  research  can  be  applied  to  accurately  characterize  the  time-frequency  characteristics  of  sea 
clutter  to  enhance  the  development  of  more  effective  clutter  mitigation  techniques,  particularly  in  the  area  of 
waveform  design  [16,  17].  This  approach  attempts  to  capture  enough  of  the  essential  signal  characteristics 
to  be  useful  while  retaining  the  low  computational  burden  of  the  statistical  models  from  which  it  is  built  on. 

1.2  Literature  Survey 

Like  many  early  developments  in  RF/microwave  technology,  the  first  systematic  presentation  of  research 
on  sea  clutter  was  published  in  a  MIT  Radiation  Laboratory  series  book  chapter  by  Ken-  and  Goldstein  [18]. 
In  this  work,  the  concept  of  a  normalized  radar  cross-section  (NRCS),  <7°,  is  introduced  along  with  a  rather 
comprehensive  description  of  sea  clutter  generation  mechanisms  as  they  were  understood  at  the  time.  Other 
early  research  discusses  the  dependence  of  the  NRCS  on  frequency,  grazing  angle,  and  wind  speeds  [19,  20]. 
Other  early  work  studied  the  scattering  from  rough  and  periodic  structures  on  the  sea  surface  [21,  22].  A 
good  summary  overview  of  the  work  of  these  early  decades  is  given  by  Barton  [23].  Much  work  focuses  on 
the  unique  characteristics  of  clutter  as  seen  from  low  grazing  angles  [24] . 

Sea  clutter  modeling  techniques  generally  fall  into  one  of  two  main  categories-  statistical  approaches 
and  quasi-deterministic  electromagnetic  approaches.  (Granted,  most  methods  arc  not  exclusively  one  or  the 
other,  but  rather,  we  arc  assigning  them  to  the  group  with  which  their  approach  most  closely  aligns.) 

Statistical  models  model  the  sea  clutter  return  as  a  stochastic  process  whose  magnitude  return  follows 
a  specified  distribution.  Rayleigh  clutter  is  considered  the  simplest  case,  arising  from  central  limit  theorem 
arguments  when  the  number  of  scatterers  in  a  resolution  cell  is  sufficiently  large  [25].  However,  when  the 
radar  resolution  increases  and  particularly  when  the  grazing  angle  is  less  than  10°  (which  is  a  common  def¬ 
inition  for  “low  grazing  angle”),  individual  wave  “spikes”  become  resolvable  and  the  distribution  becomes 
more  heavy-tailed.  While  a  number  of  distributions  have  been  proposed  to  model  low-grazing  angle  clut¬ 
ter,  one  of  the  most  prominent  was  the  K-distribution,  first  proposed  by  Jakeman  and  Pusey  [26-29]  and 
promoted  by  Ward,  Tough,  and  Watts  [25,  30-32].  Statistical  methods  have  the  advantage  that  they  arc 
relatively  computationally  inexpensive. 

Quasi-deterministic  approaches  generally  begin  by  creating  a  random  sea  surface  generated  from  a  spec¬ 
ified  wave  spectrum,  such  as  is  shown  in  Figure  1  from  which  a  deterministic  electromagnetic  scattering 
computation  is  performed  [33].  Ott  gives  a  good  overview  and  comparison  of  many  such  models,  including 
the  well-known  TEMPER  and  VTRPE  models  [34] .  While  these  methods  may  more  accurately  model  the 
fundamental  physics  of  the  problem,  they  generally  impose  a  much  more  severe  computational  burden  on 
any  sea  clutter  simulation,  which  is  a  severe  disadvantage. 

2.  MATHEMATICAL  BACKGROUND 

In  this  section  we  will  derive  statistical  characterizations  of  the  impulse  response  of  the  sea  clutter 
scattering  channel  by  modeling  it  as  a  zero-mean  WSSUS  random  process. 
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Fig.  1:  Example  swell  wave  process  realization  generated  from  a  narrow-band  gravity  wave  spectrum  from 
[33], 


airborne  radar  at  origin 


surface 


Fig.  2:  Flat  earth  geometry  in  the  x  —  z  plane. 
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(c)  Transmitted  signal  to  the  next  isorange  ring,  at  ground 
range  g  +  dg. 


(d)  Reflected  signal  from  the  next  isorange  ring,  at 
ground  range  g  +  dg. 


Fig.  3:  Illustration  of  sequential  returns  from  isorange  rings  on  the  sea  surface. 


2.1  Simulation  Geometry 

First  we  must  define  the  simulation  geometry.  We  will  assume  that  the  radar  of  interest  being  modeled 
is  in  an  airborne  platform  with  altitude  h  above  the  sea  surface,  as  shown  in  Figure  2.  In  this  figure  the 
variables  r,  g,  and  h  represent  slant  range,  ground  range,  and  altitude,  respectively,  and  9  and  a  represent 
the  depression  and  grazing  angles,  respectively.  For  a  flat  earth  geometry,  a  =  6.  The  coordinate  axes  arc 
oriented  such  that  the  platform  is  located  at  the  origin  and  the  x  —  z  plane  is  oriented  such  that  the  platform 
velocity  vector  v  =  [vx,vy,  vz]T  always  lies  in  the  x  —  z  plane.  (This  implies  that  vy  =  0.)  Because  we  are 
seeking  to  model  the  channel  characteristics  as  a  function  of  t  and  r,  we  will  express  all  geometric  quantities 
in  terms  of  these  variables.  (The  picture  in  Figure  2  is  given  for  a  fixed  value  of  t—  we  will  account  for  the 
platform  velocity  at  a  later  point.)  Thus,  the  range  r  is  given  as  r  =  \ct  and  g  is  solved  for  using  the 
Pythagorean  theorem: 


where  c  is  the  speed  of  light  and  h  is  assumed  fixed. 

We  will  model  the  radar  signal  as  a  spherical  wave  emanating  from  the  source  located  at  the  origin  as 
is  shown  in  Figure  3.  Each  small  segment  of  the  wave  reflects  successively  off  isorange  rings  on  the  sea 
surface;  the  total  return  is  therefore  the  summation  of  the  returns  from  each  isorange  ring. 
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The  power  gain  as  a  function  of  delay  can  be  obtained  using  the  radar  range  equation  [35]: 

Pr  _  G2(cj>(T),d(T ))  A2  ct0(q(t))  dA{r) 

Pt  (47r)3(^cr)4 


where 


•  Pj  =  transmitted  power  (average  power  during  pulse  duration) 

•  Pr  =  receive  power 

•  cf),  6  =  azimuth  and  depression  angles,  respectively,  such  that: 

x  =  r  cos  6  cos 
y  =  r  cos  0  sin  f 
z  =  rsin# 

•  a  =  grazing  angle,  which  for  a  flat  earth  model  is  a  =  6  =  sin_1(/i/(^cr)) 

•  G((j),  6)  =  antenna  power  gain  as  a  function  of  look  direction 

•  A  =  wavelength  =  c/f 

•  <7°  (cc)  =  normalized  RCS  (NRCS)  of  sea  surface  as  a  function  of  grazing  angle  a  for  some  fixed 
wavelength  and  sea  state 

•  dA(r)  =  differential  surface  area  of  isorange  ring  at  delay  r. 


The  differential  surface  area  dA  computation  (assuming  a  flat  earth  model)  is  illustrated  in  Figure  4  and 
is  derived  as  follows: 


dA  =  gdgd(j)=  (g)  ■  (  ^  dr  )  •  (27 r)  = 


(2tt) 


(5) 

(6) 


2.2  Wide-Sense  Stationary/Uncorrelated  Scattering  (WSSUS)  Linear  Time-Varying  (LTV)  System 
Theory 

2.2.1  Overview  and  fundamental  definitions 

The  theory  of  random  WSSUS  LTV  systems  was  first  illustrated  by  Bello  in  his  landmark  paper  [2] 
and  his  approach  has  been  used  successfully  for  decades  in  the  wireless  communication  community  to 
model  signal  amplitude  fluctuations  in  mobile  radio.  (The  mobility  of  the  radio  nodes  creates  changing 
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geometry  of  the  wireless  channel  and  makes  the  channel  impulse  response  appeal-  random.)  Because  an 
LTV  system  implicitly  models  motion,  LTV  system  theory  is  well-suited  to  modeling  the  Doppler  effect  of 
moving  scatterers,  such  as  is  present  in  a  radar  system.  In  fact,  the  only  substantial  differences  between 
the  communications  application  versus  the  radar  application  from  a  linear  systems  perspective  is  that  in 
a  radar  the  transmitter  and  receiver  are  co-located,  and  in  a  radar  system  the  antenna  is  generally  highly 
directive,  whereas  in  mobile  radio  the  antenna  is  usually  assumed  to  be  omnidirectional  [1],  (In  the  case 
of  a  single-antenna  system  this  is  true;  newer  multiple  input-multiple  output  (MIMO)  radio  systems  with 
multiple  antennas  at  both  the  base  station  and  mobile  node  violate  this  assumption.) 

To  apply  the  WSSUS  assumption  to  a  radar  scattering  channel,  we  will  express  the  impulse  response  in 
the  following  form: 


h(r,  t) 


t](t,  p )  e^2npt  dp 


(V) 


where  p  is  the  Doppler  shift  and  //(r,  p)  is  known  as  the  “delay-doppler  spread  function”  [10,  15].  This 
function  describes  the  scattering  amplitude  at  a  given  delay  (range)  and  Doppler  shift.  Substituting  (7)  into 
(2)  yields: 


y(t) 


t](t,  p)x(t 


t )  e i27Tpt  dp  dr, 


oo  J  —  oo 


(8) 
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where  x(t)  is  the  transmitted  signal  (at  baseband)  and  y(t)  is  the  received  signal.  In  (8)  it  is  clear  to  see  that 
y(t)  can  be  expressed  as  the  superposition  of  delayed  (by  r),  frequency-shifted  (by  p),  and  scaled  (by  the 
complex  gain  t](t ,  p))  copies  of  x(t).  It  should  be  noted  that  this  definition  is  general  enough  to  include  all 
scatterers  in  the  radar’s  field  of  view-  clutter  as  well  as  useful  targets.  In  this  case  the  overall  delay-Doppler 
spreading  function  is: 


^/total  —  ^targets  T  flclutter 


(9) 


due  to  the  linearity  of  the  system. 

In  the  case  of  sea  clutter,  the  multitude  of  minute  variations  of  the  sea  surface  make  it  appear  extremely 
random,  and  thus  we  will  assume  that  rj  and  h  arc  in  general  a  stochastic  process  in  this  work.  To  characterize 
the  behavior  of  sea  clutter  we  will  assume  that  it  is  zero-mean  and  compute  its  second-order  statistics: 

Rh(r,  t',  t,  t  +  At)  =  E  [h*(r,  t )  h(r' ,  t  +  At)]  (10) 

Here  we  apply  two  fundamental  assumptions:  1)  The  signal  is  wide-sense  stationary  (WSS),  implying  that 
the  autocorrelation  depends  only  on  At  on  the  time  axis,  and  2)  the  scattering  at  different  lags  r  is  uncorre¬ 
lated: 


E  [h*(r,  t )  Ii(t',  t  +  At)]  =  Rh{r,  t',  At) 

=  Ah(r,  At)  S(t  -  t')  (11) 

The  function  Ah(r,  At)  is  known  as  the  channel  correlation  function,  and  describes  the  distribution  of  the 
radar  return  as  a  function  of  lag  t  and  correlation  time  At  [1,3], 

In  an  airborne  pulse-Doppler  radar  system  we  can  think  of  the  channel  as  being  roughly  constant  during 
a  single  pulse  duration,  but  that  the  channel  changes  slightly  between  pulses  (e.g.  a  small  phase  rotation  due 
to  the  Doppler  offset  of  the  return).  Under  this  view,  we  can  treat  the  variable  r  as  the  downrange  delay  (or 
the  “fast  time”  in  a  radar  data  cube  [17])  and  the  variable  At  as  capturing  the  evolution  of  the  channel  over 
“slow  time”.  By  taking  the  Fourier  transform  of  ,4/f  (r,  At)  with  respect  to  At,  we  can  view  the  Doppler 
spectrum  of  the  return  as  a  function  of  delay.  The  so-called  “scattering  function”  is  defined  as  follows: 


/OO 

Ah(r,  At)e~j2npAtdAt. 

-OO 


(12) 


The  scattering  function  provides  a  measure  of  the  output  channel  power  as  a  function  of  delay  and  Doppler. 
It  was  first  introduced  formally  by  Price  and  Green  [5,  6]  but  was  first  hinted  at  in  a  paper  by  Westerfield 
[36].  This  function  played  a  prominent  role  in  Bello’s  fundamental  paper  where  he  coined  the  WSSUS 
nomenclature  and  showed  the  two-dimensional  Fourier  transform  relationships  (i.e.  along  the  At  yy  p 
dimension  and  along  the  r  B-  A /  dimension)  between  four  related  functions  that  equally  capture  the 
time-frequency  characteristics  of  a  WSSUS  channel,  two  of  which  we  have  described  here  (,4/,(r,  At)  and 
Cn ( t,  p)).  As  was  mentioned  earlier,  the  scattering  function  has  been  used  to  characterize  radar  returns  since 
its  formal  inception  by  Price  and  Green,  who  used  it  to  model  returns  from  large  distributed  targets  such  as 
the  return  from  a  planet  in  radar  astronomy.  Van  Trees  also  briefly  describes  the  scattering  function  in  his 
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most  well-known  book  [7],  as  does  Kay  in  at  least  one  paper  [8],  but  despite  the  prestige  associated  with 
some  of  its  proponents,  it  has  remained  a  niche  mathematical  tool  that  sees  little  practical  use.  It  is  uncertain 
why  this  is  the  case,  given  the  relative  simplicity  it  allows  one  to  describe  a  radar  scattering  channel. 

It  can  also  be  shown  that  the  autocorrelation  function  of  the  delay-doppler  spread  function  rj  for  a  WS- 
SUS  channel  is: 


®  [v*(,Ti  p)  v(T'i  P)\  =  Cv(t,  p)  S(t  —  t')  S(p  —  p').  (13) 

(It  is  sometimes  easier  to  work  with  this  form  when  deriving  new  expressions.) 

Note  that  nowhere  in  this  discussion  have  we  yet  specified  a  distribution  for  the  samples  of  1i(t.  t).  If  we 
take  the  further  step  of  specifying  that  h  is  a  complex  Gaussian  random  process  with  independent  real  and 
imaginary  components,  then  the  distribution  of  h  is  completely  specified  by  either  the  correlation  function 
Ah  or  the  scattering  function  Cv,  and  sample  realizations  of  h  can  be  generated  by  filtering  white  noise 
to  have  the  Doppler  spectrum  determined  by  Cv.  It  is  well  known  that  when  the  values  of  h  are  given  by 
a  complex  Gaussian,  the  magnitude  envelope  \h\  will  be  Rayleigh-distributed  and  the  power  \h\2  will  be 
exponentially  distributed  [1].  The  Rayleigh  distribution  was  the  de-facto  standard  for  modeling  sea  clutter 
returns  for  early  radars  and  is  a  good  approximation  when  the  beam  footprint  (and  range  resolution)  is  wide 
and/or  the  grazing  angle  is  large  and  thus  encompasses  a  large  number  of  independent  scatterers  whose 
returns  are  assumed  to  sum  to  be  approximately  Gaussian  using  central  limit  theorem  arguments  [25,  37]. 

2.2.2  Validity  of  the  WSSUS  assumption 

We  have  assumed  that  the  channel  is  both  wide-sense  stationary  (WSS)  (statistics  depend  only  on  At, 
not  t)  and  that  there  is  uncorrelated  scattering  (US)  (i.e.  the  returns  are  uncorrelated  for  different  values 
of  r).  The  WSS  assumption  is  relatively  noncontroversial-  over  the  course  of  a  single  CPI,  which  is  on 
the  order  of  milliseconds,  or  even  several  CPIs,  it  seems  very  reasonable  that  the  physical  state  of  the  sea 
will  have  not  changed  significantly.  (This  is  equivalent  to  saying  that  the  joint  probability  density  of  the  sea 
surface  properties  retain  the  same  statistics  over  this  time  period.)  The  US  assumption,  however,  is  much 
more  controversial,  due  to  the  fact  that  the  dimensions  of  larger  sea  surface  features  (such  as  swell  waves) 
are  much  larger  than  the  radar  resolution,  implying  that  there  will  be  correlation  between  the  returns  in 
closely-spaced  range  bins  [25]. 

One  way  of  introducing  correlation  in  the  r  dimension  would  be  to  replace  the  delta  function  in  (11) 
with  a  function  such  as  a  decaying  exponential  to  express  correlation  between  adjacent  range  samples: 


Rh(T,Tf,  At)  =  Ah(T,  At)e  |r  r'l/Tc, 


(14) 


where  Tc  is  a  correlation  length,  usually  determined  empirically.  Many  models  exist  for  modeling  spatial 
correlation  of  sea  clutter  [38];  one  formula  for  correlation  length  (in  meters)  is  given  by  [24]: 


7 r  U2 


=  —  —  3  cos2  4>w  +  1  =  ~cTc 

2  g  V  Y  2 


(15) 
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where  U  is  the  wind  speed  in  m/s,  g  ~  9.81  m/s  is  the  gravitational  constant,  and  <pw  is  the  wind  direction. 
For  example,  if  the  wind  speed  is  10  m/s  (corresponding  approximately  to  Douglas  sea  state  4  or  “rough 
sea”  [25]),  the  maximum  correlation  length  is  32.02  m  or  approximately  213.6  ns  of  lag.  If  this  distance 
is  greater  than  the  range  resolution  of  the  radar  operating  in  that  environment,  then  it  cannot  be  assumed 
that  the  clutter  return  is  white  and  matched  filter-based  detectors  will  perform  sub-optimally  [17].  In  other 
words,  estimates  of  the  signal-to-clutter  ratio  based  on  the  WSSUS  assumption  will  be  too  optimistic  if 
significant  range  correlation  exists.  One  way  of  modeling  the  spatial  correlation  of  clutter  would  be  to 
generate  WSSUS  clutter  samples  then  filter  them  along  the  r  axis  with  a  filter  designed  using  (15)  whose 
spectral  characteristics  approximate  that  of  the  clutter  spectrum. 

Non-WSSUS  LTV  systems  have  been  studied  by  some  authors  (e.g.  [39]),  but  for  simplicity’s  sake  we 
will  use  the  WSSUS  model  in  this  work  so  we  can  produce  an  elementary  result  to  fill  the  immediate  need  for 
a  first-order  model  to  develop  further  intuition  of  sea  clutter  scattering  channels  and  leave  the  non-WSSUS 
formulation  to  future  work. 


3,  WSSUS  SEA  CLUTTER  SCATTERING  FUNCTION 

To  create  a  model  for  Cv(t.  p),  we  need  to  characterize  the  signal  return  in  range  (via  r)  and  Doppler 
(via  p).  The  majority  of  the  decay  in  signal  amplitude  with  r  we  have  already  accounted  for  in  (4).  To 
account  for  the  Doppler  characteristics  we  will  consider  the  return  from  a  single  isorange  ring  and  then  later 
combine  our  results  into  a  single  function. 

3.1  Delay-Doppler  Derivation 

Using  the  Clarke  model  for  uniform  scattering  [40, 41]  that  is  widely  used  in  the  wireless  communication 
field  as  a  guide,  we  will  derive  Doppler  spectrum  of  the  scattered  sea  clutter  signal  at  a  single  small  isorange 
ring  located  at  ground  range  g  with  width  dg. 

We  will  assume  that  the  impulse  response  from  this  ring  is  a  sum  of  N  delta  functions  with  a  time- 
varying  amplitude  and  phase  (which  encompasses  the  Doppler  shift  as  well).  We  will  segment  the  domain 
of  the  ring  into  evenly-spaced  angular  patches  in  azimuth  d  £  [0,  2n)  with  random  amplitudes  and  phases 
as  shown  in  Figure  5.  We  will  also  assume  that  the  scatterers  on  this  ring  are  located  at  delay  r'  and  write 
the  impulse  response  as  follows: 


N- 1 

h(r,t',  t')  =  ^2  an(t)e~jxl’ri^S(T  -  r') 

n=0 


(16) 


where  an  is  the  real- valued  amplitude  of  scatterer  n  and  the  phase  of  scatterer  n  is  given  by: 

=  27 r/r  -  27 xpnt  -  (n,  (17) 

where  /  is  the  carrier  frequency,  r  and  pn  are  the  delay  and  Doppler  shift,  respectively,  of  scatterer  n,  and 
Cn  is  a  uniformly-distributed  random  phase  over  [0,  2n).  We  will  assume  that  the  amplitude  of  each  return 
does  not  change  significantly  over  a  coherent  processing  interval  (CPI),  therefore  an(t)  ~  an.  We  will  also 
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Fig.  5:  Clutter  patch  returns  for  a  single  isorange  ring. 


assume  that  an  and  arc  independent  of  each  other.  Note  that  because  the  delay  t'  to  each  scatterer  is  the 
same  and  does  not  depend  on  n  (which  represents  the  azimuth  angle),  we  can  factor  it  out  of  the  sum: 


JV— t 

h(r,  t;  t')  =  5(t  —  t')  Y^  ane~^n^ 

n= 0 


(18) 


The  autocorrelation  function  of  the  clutter  signal  is  given  by: 


Rh{n,T2,t,t  +  At)  =  E  [h*(Ti,t)  h(r2,t  +  At)] 


=  E 


N—l 


N—l 


$(ti-t')  Y 


&r)  €■' 


5(r2  -  t')  Y 


nme 


n= 0 


m= 0 


N—l  N—l 


S(ti  -  t')8(t2  -  t')  EE1 


m= 0  n= 0 


QjmQ,n£ 


j (ipm  (t+At)-lj>n  (t)) 


(19) 
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We  will  assume  that  the  phase  shifts  {f>n}  arc  mutually  independent,  therefore: 


TV— 1  TV-1 

Rh{n,T2,t,t  +  At)  =  S(t i  -  t')5(t2  -  t')  EE  E  [aman]  E 


m= 0 n= 0 


N—l 


S(ti  -  t')5{t2  -  t')  ^2  E  [all  e^W"(t+A*)_V,"(t)) 


n= 0 


TV-1 


<S(ri  —  t')5(t2  —  t')  ^  E  [al]  e'2'KPr 


At 


n= 0 


=  Rh{ri,  t2i  At]  t') 


(20) 


Note  that  E  [al]  is  simply  the  channel  gain  from  scatterer  n,  therefore  using  (4)  and  scaling  it  to  account 
for  the  fact  that  the  surface  area  is  smaller  by  a  factor  of  1  /N,  we  obtain: 


ttt<  r  2i  _  1  Rr 

E  ^  ~  NPt 


G2(nA<fi,9)  A2  </*_  _  7T  2  ,  ,  , 
(47t)3  (2ct')4  N  2 


(21) 


Note  that  N  =  2tt/A6,  where  A <4  is  the  angular  spacing  between  patches,  which  upon  substituting  in  (21) 
yields: 


A2<r°  •  f  c2t' dr’ 
(4vr)3  (^ct')4 


—G2(nA(j),  9)  A<f 
27 r 


(22) 


Combining  this  with  (20)  and  taking  the  limit  as  N  — >  oo  yields 


,  ,  ,  A2o-ufc2r  or' 

t-2,  At;  r  )  =  d(n  -  r  )<5(t2  -  t  ) 


TV-1 


=  5(n  -  t’)5{t2  -  r') 


i - ,  lim  —  G2{nAcj),  0)  < 

(47 r)3  (jcr')4  TV-roo  27r 

aAV°  •  fcVdr'  1  /-  ^  m  j2,P’^e)At 


J2irpnAt 


Af 


(47r)3  (her')4  27 r  J_- 


G2(4>,  9)  e3 


dej) 


(23) 


The  Doppler  shift  of  the  infinitesimal  scatterer  located  at  azimuth  cj>,  p'(4>),  is  given  as  follows: 


/  ^ 

p  (0,  9)  =  —  (yx  cos  9  cos  f>  +  vy  cos  9  sin  <f>  +  vz  sin  9) 
A 

=  px  cos  9  cos  4>  +  py  cos  9  sin  +  pz  sin  9 


(24) 


where  vx,  vy,  and  vz  are  the  x,  y,  and  z  components  of  the  aircraft  velocity  vector,  respectively,  and 
px  =  2vx/X ,  py  =  2vy/X,  and  pz  =  2vz/X.  Note  that  this  autocorrelation  function  corresponds  to 
the  return  from  a  single  isorange  ring;  to  obtain  the  autocorrelation  from  the  entire  sea  surface  we  will  take 
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advantage  of  the  US  assumption  of  this  channel  that  the  return  from  different  range  increments  arc  uncorre¬ 
lated  and  thus  that  the  overall  autocorrelation  is  just  the  sum  of  the  autocorrelations  from  each  ring,  implying 
an  integral  over  t'\ 


Rh{Ti,T2,At )  =  j  Rh{r\,T2,  At;  t'). 


(25) 


Note  the  omission  of  the  differential  dr'  due  to  it  already  being  included  in  Rh(Ti,T2,  At;  t').  Using  the 
sifting  property  of  the  delta  function  to  complete  the  integral  yields: 


Rh{Ti,T2,  At;  t')  =  5[t i  -  r2) 


A2c7°  •  f  c2n  l 


(47r)3  (5CT1)4  27 t  J_ 

=  5 (n  -t2)  Ah(n,  At) 


G2(cj>,9)e 


j27rpf  ((f), 0)  At 


dcj) 


(26) 


Note  that  the  autocorrelation  function  is  in  the  form  required  for  a  WSSUS  random  process  and  thus  we  can 
immediately  identify  the  WSSUS  correlation  function  ,4/,(r.  At)  from  (26): 


a-{t’  At) = ^  L s)  ei2rpW)A‘  "* 


A  2C(T°  1 


(47r)3(2cr)3  2 


G2((j),d)ej2np'^e)At  dct> 


A  2ca° 


(47 r)3  {\ct )3  Jo 


G2(0,  0)ej2*S(M)At 


(27) 


where  the  last  equality  is  due  to  the  fact  that  the  integrand  is  an  even  function  of  (f>.  To  find  the  scattering 
function  Cv,  we  take  the  Fourier  transform  of  (27)  with  respect  to  At: 


Cv(t,  p)  =  J^At  [■ Ah(r ,  At)] 


\2ca° 


(47r)3  (^cr)3  Jo 
A  2ca° 


G2^0)FAt 


oj2np'(4>,0)At 


dcj) 


(47r)3  (icr)3  Jo 


G2((t>,  0)  5{p  -  p((j),  9))  dcj) 


(28) 


The  integral  in  (28)  must  be  solved  using  a  change  of  variable  to  eliminate  the  0  dependence  in  the  argument 
of  the  delta  function  to  allow  for  the  evaluation  of  the  integral.  Recall  that  under  our  coordinate  system 
definition  vy  =  0,  therefore  the  Doppler  shift  is: 


p(<t>,  0)  =  px  cos  9  cos  4>  +  pz  sin  9 
p>  -  pz  sin  9 


(29) 


COS  (j)  = 


px  cos  9 


(30) 


Cn{r,p)  —  Cv(r,p) 
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ig.  6:  Normalized  Doppler  spectrum  for  an  isotropic  radiator  with  vz  =  0. 


Fig.  7:  Normalized  Doppler  spectrum  for  a  directive  radiator  with  vz  =  0. 
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Over  the  region  of  integration,  we  can  solve  for  0  and  perform  the  variable  substitution: 


dtp  =  — 


1 


=  cos 


1 


_1  p  -  Pz  sm( 


Px  COS  6  1  /  ,  -0X2 

1  I  P  ~PZ  sln^ 

px  cos  9 


px  COS  0 

dp  = 


-1 


y f {px  cos  9)2  -  {p'  -  pz  sill  6) 


-.dp  (31) 


G2(4>,0)  5(p  -  p((j),9))d<p 


f* pz  sin  6-\-px  cos  6 


g2(4>(p'),  0)  5  [P  —  P') 


t pz  sin  9  px  cos  9  \/  (pxCOS0)2  —  (p'  —  pzSUl8)2 


dp1 


=  < 


G2(<Kp),9) 


\/ ( px  cos  9)2-(p-pz  sin 9)2  ’ 

0,  else. 


|  p  —  pz  sin  6 1  <  px  cos  6 


(32) 


Combining  this  result  with  (28)  yields: 


Cv(t,p)  =  < 


A2 


C  (7 


(4tt)3  (jCt)3  yj cos0)2_(p_p,  sin 9) 


0, 


G~J?{p)'9)  ■  2 ?  \p- pzsin9\  <  Pxcost 


else. 


(33) 


Many  of  the  factors  in  (33)  have  an  implicit  dependence  on  r  that  has  been  suppressed  for  brevity-  the 
full  form  is  shown  below: 


Cv(t,  p) 


(  A 2c  .  <7°(a{r))  _  _ G2((j>(T,p),9(r)) _ 

(iiTp"  (|cr)3  y/ cos(6(T)))2-(p-pz  sin (0(r)))2  ’ 


|p  -  p2sin(0(r))|  <  pa;Cos(6l(r)) 

(34) 


0, 


else. 


This  result  generalizes  the  Clarke  model  to  the  case  when  the  scatterers  are  no  longer  located  in  the  same 
plane  as  the  velocity  vector  of  the  radar  [40].  A  plot  of  the  normalized  Doppler  spectrum  for  a  given  range 
was  generated  using  (32)  and  is  plotted  in  Figure  6  for  an  isotropic  radiator  (i.e.  G  =  1)  and  for  a  directive 
radiator  in  Figure  7  with  antenna  pattern  G(4>,  9)  =  3sinc2(5(</>  —  f  ))  (i.e.  the  antenna  is  pointed  45° 
counterclockwise  in  azimuth  from  the  a;— axis).  It  can  be  seen  from  Figure  6  that  when  the  energy  is  evenly 
radiated  in  all  directions  the  Doppler  shift  of  the  returns  is  focused  near  ± px  cos  9,  the  peak  Doppler  shift, 
but  it  can  be  seen  from  Figure  7  that  the  directivity  of  the  antenna  focuses  the  Doppler  spectrum  around  the 
Doppler  shift  of  the  patch  of  the  sea  surface  the  beam  is  pointed  at. 
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Fig.  8:  Scattering  function  of  sea  clutter  channel. 


A  plot  of  the  scattering  function  over  the  full  r  —  p  plane  was  generated  with  a  circularly- symmetric 
sine2  antenna  pattern  generated  as  follows.  The  basic  antenna  pattern  shape  (pointed  at  (j)  =  9  =  0  rad)  is 
given  by: 


+  ,  (35) 

where  D  is  the  directivity  in  linear  units,  B  is  the  beamwidth  in  rad,  W  =  0.88589  rad  is  a  normalization 
constant,  and  (d,  9)  are  given  in  rad.  Correctly  rotating  the  antenna  pattern  to  depression  angle  9a  requires 
a  nonlinear  transformation  of  the  angular  coordinates: 

G(4>,  9)  =  GoWW,  9 ■  9a),  9'(<t>,  9 ;  9a))  (36) 

with  the  transformation  equations  given  as  follows: 

<t>'(<j),  9 ;  9  a)  =  atan2(cos  9  sin  cos  9a  cos  9  cos  +  sin  9a  sin  9)  (37) 

9>(<j),  9;  9 a)  =  atan2(—  sin  9a  cos  9  cos  (j)  +  cos  9a  sin  9 , 
y/  (cos  9  a  cos  9  cos  (?)  +  sin  9a  sin  9)2  +  (cos  9  sin  <f)2)  (38) 


Go((j>,  9)  =  D  sine 


A  plot  of  the  full  scattering  function  Cv(t,  p)  is  computed  and  plotted  in  Figure  8  for  /  =  10  GHz, 
h  =  5  km,  v  =  [300, 0, 0]T  m/s,  9a  =  25°,  B  =  10°  =  0.17453  rad,  and  D  =  100  =  20  dB.  NRCS 
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(<t°)  data  for  VV  polarization  at  X-band  as  a  function  of  grazing  angle  a  was  obtained  by  interpolating  data 
points  taken  from  Figure  7.13  from  [37]. 

3.2  Time-Frequency  Power  Distribution  at  Channel  Output 

One  of  the  goals  of  pulse-Doppler  processing  is  to  use  Doppler  filtering  to  isolate  a  target  in  delay- 
Doppler  space  and  obtain  improved  detection  performance  [35].  We  can  apply  the  theory  we  have  developed 
for  WSSUS  scattering  channels  and  use  the  scattering  function  Cv  to  predict  the  time-frequency  power 
distribution  of  the  output  signal.  This  presentation  is  essentially  an  updated  version  of  the  derivation  given 
in  [6]. 

We  will  first  use  the  definition  of  the  radar  ambiguity  function  from  [42]  as  our  building  block: 


/OO 

u{t)u*(t  +  T)ej2nptdt., 

-OO 


(39) 


which  can  be  thought  of  as  the  matched  filter  output  of  a  unit-energy  radar  pulse  u{t)  delayed  in  time  by 
r  and  Doppler-shifted  in  frequency  by  p.  Consistent  with  our  earlier  definitions,  we  will  note  that  the 
transmitted  signal  x(t )  is  expressed  in  terms  of  u(t)  as  follows: 


x{t)  =  [fi 

The  channel  output  before  matched-filtering,  sampled  at  time  r  is  given  by: 

r(r)  =  J  j  r){r' ,  p')x(r  —  r/)e?27rp  T cLt' dp1 . 


\x{t')\2dt'  J  u{t)  =  \/Rxu(t) 


(40) 


(41) 


If  we  consider  just  the  returns  from  a  single  Doppler  shift  p  and  downshift  them  to  baseband  for  observation 
we  obtain: 


?’(r;  p)  =  e  •?27rpT  J  J  t](t' ,  p')x(r  —  r/)e?27rp  Tdr' dp' 
ri(r\  p')x(t  —  r/)e-?27r0’  ~P^T dr' dp' . 


(42) 


This  can  be  thought  of  as  the  output  from  a  Doppler  filterbank  filter  tuned  to  p.  The  matched-filtered  output 
is  given  by: 


y(r-,p)  =  u*(—r)  *  r(r;  p) 

=  J  j  j  t](t' ,  p')x{t'  —  T/)ej27r^p  ~pS>t  u*(t'  —  r)dt' dr' dp' 

=  \Z~Ex,  f  f  f  y(T' )  p')u{t'  —  T,)ei2'K(p  ~p^  u*(t'  —  r)dt' dT1 dp' . 


(43) 
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Fig.  9:  Ambiguity  function  for  a  train  of  unmodulated  rectangular  pulses. 


Let  t!  —  t'  =  t  and  substitute  this  into  (43),  which  yields: 


y{r;p)  =  y/E. 


=  \E, 


/»«>*(*  +  r'-ry^mrdp' 

(r,  p)ej27T(p'~p)T'  x{t  —  t,  p  —  p)dr  dp  . 


(44) 


The  filter  output  of  (44)  is  a  random  process  due  to  77;  to  characterize  the  output  power  distribution  we 
will  take  the  expectation  of  the  magnitude  squared: 


=  Er 


P{t,p)  =  E  [\y(r;  p)|2] 

E  e-V*Wi-p)<ei2*te-PK 

-  T,  p'l  -  p)x{r2  -  T,  f>2  -  p)dT[dT2dp1dp2. 


(45) 


Using  (13),  we  obtain: 


P(r,p)  =  Ex  J  J  Cv(t' ,  p')\x{t'  -r,p'  -  p^dr'dp 

=  Ex  Cv (r ,  p)  *  *\x(~T,-p)\2 
=  Ex  Cv(t,  p)  *  *|x(t,  p) I2 


(46) 
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Fig.  10:  Output  power  distribution  from  a  sea  clutter  channel  with  an  input  pulse  train  of  unmodulated 
rectangular  pulses. 


where  **  denotes  two-dimensional  convolution,  and  the  identity  | %(— r,  —p)\  =  \x(T>  p)\  has  been  used  to 
simplify  the  expression. 


To  illustrate  an  example  scattered  power  distribution,  we  will  use  a  train  of  constant-frequency  rectan¬ 
gular  pulses  of  duration  T.  The  equation  for  a  single  pulse  is  given  by: 


u(t) 


(47) 


whose  ambiguity  function  magnitude  squared  is  given  by  [42] : 


I Xt(t,  p) I2 


(l  -  sine  (Tp  (l  -  ^)) 


t\  <  T 


(48) 


0 


else. 
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The  ambiguity  function  magnitude-squared  of  the  train  of  N  pulses  with  pulse  repetition  interval  (PRI)  Tr 
is  given  by: 


\x(r,p)\2  = 


Tf  £p=-(JV-l)  \Xt(t  -pTr,p)\ 


sin(7rp(Ar— |p|)Tr) 


sin(npTr) 


<  NTr 


(49) 


0, 


else. 


The  magnitude-squared  ambiguity  function  of  the  pulse  train  is  plotted  in  Figure  9  for  a  train  with  N  =  4 
pulses  with  pulse  width  T  =  10  ps  and  PRI  Tr  =  40  ps. 


Using  equation  (46)  and  the  scattering  function  plotted  in  Figure  8  we  compute  the  output  power  dis¬ 
tribution  P(t,  p)  and  plot  it  in  Figure  10.  It  can  be  seen  that  the  output  power  distribution  consists  of  a 
repetition  of  blurred  copies  of  the  scattering  function  Crj. 


3.3  Signal-to-Clutter  Ratio 


Because  the  output  of  the  radar  scattering  channel  is  linear,  we  can  use  use  the  superposition  principle 
of  (9)  to  model  the  return  from  multiple  objects,  such  as  a  ship  in  sea  clutter.  The  scattering  function  for  a 
point  target  with  RCS  cr  at  delay  t  and  Doppler  shift  p  is  derived  using  the  radar  range  equation: 


Cv'{t,  p) 


(47 r)3(4cr)4 


6(t 


t)5(p-p), 


(50) 


which  implies  that  the  output  power  distribution  is  given  by: 


P'(r,p)  =  ExCv>(t,p)  **| x{T,p)f 


Ex - (A  \ 3 /  1  -X4 - \x{r-T,p-p)\  . 

(47t  )d(^cr)4 


(51) 


It  is  clear  from  this  result  that  the  resolution  in  delay-Doppler  is  fundamentally  limited  by  the  resolution  of 
the  radar  pulse,  as  manifested  by  the  width  of  the  ambiguity  function  along  the  r  and  p  axes. 


Detection  probability  is  usually  expressed  as  a  function  of  the  signal-to-noise-plus-clutter  ratio  or  just 
the  signal-to-clutter  ratio  if  the  detection  performance  is  clutter-limited  [17].  Using  (34),  (46),  and  (5 1),  we 
can  now  easily  compute  the  signal-to-clutter  ratio  (SCR)  as  a  function  of  delay  and  Doppler,  7 (r,  p),  for  a 


20 


Corey  D.  Cooke 


point  target  embedded  in  sea  clutter: 


7  (t,  p)  = 


P'(r,p) 

p(r,  P ) 


=  < 


G2(^(r,p),e(r))o-|x(r— r,p— p)|2/(ter)4 


ccr®  (cx(t)) 


(hCT)3\/ Cx  cos (0(t)))2-(p~pz  sin(0(i-)))2 


**lx(r>p)l2 


oo, 


Ip  -  pzsin(0(r))|  <  px  cos (0(r)) 


else. 


(52) 


Note  that  in  the  regions  of  the  r  —  p  plane  where  the  clutter  power  is  nonzero,  the  SCR  is  independent  of 
the  transmitted  signal  energy.  This  result  is  general  enough  to  account  for  straddling  loss  if  the  spectrum  is 
sampled  at  points  r  /  f  and  p  ^  p  in  a  digital  processing  system. 


4.  CONCLUSION  AND  FUTURE  WORK 

In  this  work  we  have  presented  a  new  mathematical  formulation  into  which  the  radar  sea  clutter  modeling 
problem  can  be  cast  which  allows  for  a  conceptually  easy  way  to  study  the  time -frequency  characteristics  of 
the  scattered  signal.  By  presenting  the  problem  in  this  form,  the  received  range-Doppler  map  can  be  easily 
and  quickly  computed  for  an  arbitrary  pulse  shape,  antenna  pattern,  and  NRCS  angular  distribution. 

Future  work  will  focus  on  studying  the  validity  of  the  fundamental  assumptions  undertaken,  particularly 
the  assumption  that  the  scattering  channel  can  be  modeled  accurately  by  a  WSSUS  random  process.  It  is 
hoped  that  the  compactness  of  this  mathematical  representation  will  facilitate  more  rapid  development  of 
effective  clutter  mitigation  techniques  in  the  future. 
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